%% This file is called to compute comparative static results; 
%%%% It needs "par" as parameters
%%%% It outputs results stored in variables named similar to Inv_model_4
%%%% below

disp('-------------> Model 4, idiosyncratic shocks + imperfect reallocation + perfect credit <---------')

%% Initialization
Inv_model_4                 = zeros(length(ALPHA_grid), length(THETA_grid));
Prod_model_4                = Inv_model_4;
Prod_std_model_4            = Inv_model_4;
Output_model_4              = Inv_model_4;
Hours_model_4               = Inv_model_4;
R_ratio_model_4             = Inv_model_4;
P_share_model_4             = Inv_model_4;
Aqc_model_4                 = Inv_model_4;
Aqc_P_model_4               = Inv_model_4;
Sppe_model_4                = Inv_model_4;
Sppe_P_model_4              = Inv_model_4;
Util_model_4                = Inv_model_4;
Welfare_model_4             = Inv_model_4;
Consumption_model_4         = Inv_model_4;
Sec_mkt_price_model_4       = Inv_model_4;
Z_Y_model_4                 = Inv_model_4;
Z_K_model_4                 = Inv_model_4;
debt_model_4                = Inv_model_4;
Threshold_model_4           = Inv_model_4;
perfect_reallocation        = 0;
cali_step                   = 0;

%% Comparative statics
for ii = 1 : length(THETA_grid)
    par = par_ini;
    par.THETA       = THETA_grid(ii);
    fun             = define_function_fun(par);
    
    fprintf('I am doing the case of THETA = %.3f\n\n', THETA_grid(ii))
    
    Inv_temp        = zeros(length(ALPHA_grid), 1);
    Prod_temp       = Inv_temp;
    Prod_std_temp   = Inv_temp;
    Output_temp     = Inv_temp;
    Hours_temp      = Inv_temp;
    R_ratio_temp    = Inv_temp;
    P_share_temp    = Inv_temp;
    Aqc_temp        = Inv_temp;
    Aqc_P_temp      = Inv_temp;
    Sppe_temp       = Inv_temp;
    Sppe_P_temp     = Inv_temp;
    Util_temp       = Inv_temp;
    Welfare_temp    = Inv_temp;
    Consumption_temp= Inv_temp;
    Sec_mkt_price_temp = Inv_temp;
    Z_Y_temp        = Inv_temp;
    Z_K_temp        = Inv_temp;
    debt_temp       = Inv_temp;
    Threshold_temp  = Inv_temp;
    B0              = 0.05;
    
   parfor jj = 1: length(ALPHA_grid)
        par_temp = par;
        par_temp.ALPHA  = ALPHA_grid(jj);%
        
        x       = fsolve(@(x)eqm_ss_iid_cases_2_to_5(x, par_temp, fun, perfect_reallocation, cali_step), [B0, par_temp.ALPHA_0], options);
        B       = x(1);
        ALPHA_0 = x(2);
        
        [fval, eqm]          = eqm_ss_iid_cases_2_to_5(x, par_temp, fun, perfect_reallocation, cali_step);
        eqm.IM               = average_log_normal(eqm.surplus_net, par.c_mu, par.c_sig) / par.XI_IM;  
       
        Consumption_temp(jj) = eqm.c_star + eqm.w * eqm.IM * par_temp.ZETA;
        Inv_temp(jj)         = eqm.K * par_temp.DELTA; 
        Hours_temp(jj)       = eqm.H;            
        Output_temp(jj)      = (eqm.Y + eqm.w * eqm.IM * par.ZETA) / Y_Hansen;
        Prod_temp(jj)        = eqm.K_bar_K / K_bar_K_Hansen;
        Prod_std_temp(jj)    = eqm.prod_std / prod_std_Hansen;
        R_ratio_temp(jj)     = (eqm.E_a_K + eqm.E_p_K) / (eqm.E_a_K + eqm.E_p_K + par_temp.DELTA);
        P_share_temp(jj)     = eqm.E_p_K / (eqm.E_a_K + eqm.E_p_K);
        Aqc_temp(jj)         = eqm.E_a_K  / (eqm.E_a_K + eqm.E_p_K + par_temp.DELTA);
        Aqc_P_temp(jj)       = eqm.prob_A;
        Sppe_temp(jj)        = eqm.E_p_K  / (eqm.E_a_K + eqm.E_p_K + par_temp.DELTA);
        Sppe_P_temp(jj)      = eqm.prob_P;
        Util_temp(jj)        = log(eqm.c_star) - par_temp.XI * eqm.H;
        %Welfare_temp(jj)  = eqm.c_star / c_star_Hansen * exp(par_temp.XI * (H_Hansen - eqm.H)) - 1;
        Welfare_temp(jj)     = eqm.c_star / eqm_calibrated.c_star * exp(par_temp.XI * (eqm_calibrated.H - eqm.H) + par_temp.XI_IM * (eqm_calibrated.IM - eqm.IM)) - 1; % use calibrated c_star and H

        Sec_mkt_price_temp(jj) = eqm.average_2nd_price;
        Z_Y_temp(jj)         = ALPHA_0 * eqm.Z / eqm.Y;
        Z_K_temp(jj)         = ALPHA_0 * eqm.Z / eqm.K;
        debt_temp(jj)        = ALPHA_0 * par_temp.ALPHA * eqm.debt / eqm.Y;
        Threshold_temp(jj)   = eqm.threshold_iota;        
 
    end
    Inv_model_4(:,ii)      = Inv_temp;
    Prod_model_4(:,ii)     = Prod_temp;
    Prod_std_model_4(:,ii) = Prod_std_temp;
    Output_model_4(:,ii)   = Output_temp;
    Hours_model_4(:,ii)    = Hours_temp;
    R_ratio_model_4(:,ii)  = R_ratio_temp;
    P_share_model_4(:,ii)  = P_share_temp;
    Aqc_model_4(:,ii)      = Aqc_temp;
    Aqc_P_model_4(:,ii)    = Aqc_P_temp;
    Sppe_model_4(:,ii)     = Sppe_temp;
    Sppe_P_model_4(:,ii)   = Sppe_P_temp;
    Util_model_4(:,ii)     = Util_temp;
    Welfare_model_4(:,ii)  = Welfare_temp;
    Z_Y_model_4(:,ii)      = Z_Y_temp;
    Z_K_model_4(:,ii)      = Z_K_temp;
    debt_model_4(:,ii)     = debt_temp;
    Threshold_model_4(:,ii)= Threshold_temp;    
    Consumption_model_4(:,ii)   = Consumption_temp;
    Sec_mkt_price_model_4(:,ii) = Sec_mkt_price_temp;   
end
